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Normal states of the attractive Hubbard model, especially in two dimension, are studied 
in the light of a transition from a Fermi liquid to an insulating or gapped state. A series of 
variational Monte Carlo calculations with better statistics is carried out to estimate accu- 
rately expectation values by several many-body wave functions. Although a relatively clear 
crossover is observed even in the plain Gutzwiller wave function, the states in both regimes 
are metallic. Meanwhile, a substantial metal-insulator transition takes place at \U\ ~ W 
(band width) in an improved wave function in which intersite correlation is introduced by 
taking account of virtual states in the second-order perturbation in the infinite- \U\ limit. The 
critical value is favorably compared with recent results of the dynamical-mean-field approxi- 
mation. In contrast, a conventional Jastrow-type wave functions scarcely improve the normal 
state. In addition, the issue of Brinkman-Rice metal-insulator transition is reconsidered with 
much larger systems. 

§1. Introduction 

To consider the pseudogap phenomena found in cuprate superconductors, knowl- 
edge as to the attractive Hubbard model (AHM) is useful especially in the following 
points: i) AHM undergoes a crossover of the mechanism of superconductivity from a 
BCS type to a Bose-condensation one with increasing the correlation strength. ^ ii) 
Controllable and quantitative calculations are relatively easy, iii) Through a canon- 
ical transformation AHM on typical bipartite lattices connects with the repulsive 
Hubbard model (RHM) — a plausible model for the high-T c cuprates in two dimen- 
sion (2D). To understand the item i) as well as pseudogap, it is indispensable to study 
normal-state properties in low temperature, even though an s-wave superconducting 
phase prevails in the whole ground-state phase diagram of AHM in higher lattice di- 
mensions. 2 )' 3 ) So far, the above crossover and accompanying anomalous or gapped 
normal-state properties, for instance, in static spin susceptibility and single-particle 
spectral function have been studied by such approaches as quantum Monte Carlo 
(QMC) methods, 4 ^' 5 ^ T-matrix approximations. 6 ^ Quite recently it was shown by 
dynamical mean-field approximations (DMFT) 7 )' 8 ) that the crossover from a Fermi 
liquid to a gapped state at finite temperature is linked to a metal-insulator transition, 
which arises at \U\ ~ W (W: band width), in low temperature. 9 ^ 

Thanks to these studies, our understanding on the normal state of AHM has 
been promoted. However, many of the above results are not conclusive, because 
each of the above approaches has weak points. T-matrix approximations are reliably 
applied only to a low-density regime of weak to intermediate correlation strength. 
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The validity of DMFT is not confirmed for realistic dimensions, particularly in two 
dimension (2D), and the information on spatial correlation is not drawn. As for 
QMC methods, it is not easy to pursue a normal state for T < T c . In addition, the 
reliability deteriorates with increasing correlation strength beyond the band width 
as well as the system size, 10 ) although the minus sign problem is absent. Thus, it is 
necessary to corroborate those findings by complementary methods. 

The purpose of this paper is to study the normal-state properties of AHM by 
applying a variational Monte Carlo (VMC) method directly to 2D systems. Merits 
of the VMC approach lie in the applicability to any correlation strength and any 
electron density, to a variety of many-body trial wave functions, and to large enough 
lattices to check system-size dependence. Taking advantage of these merit, we espe- 
cially lay emphasis on crossovers or metal-insulator transitions of the states with the 
above situation in mind. As another aim we consider the resemblance and difference 
between AHM and RHM, which are interconnected by a canonical transformation. 
Such a check is indispensable before a direct comparison of the results of AHM with 
the experiments of the high-T c cuprates. 

This paper is organized as follows: In §2 the model is defined, and the rela- 
tionship between the attractive and the repulsive Hamiltonians is briefly described. 
Thereby, the previous results of RHM sometimes become applicable to AHM. In §3 
the well-known Gutzwiller wave function (GWF) is studied to show a metal-to-metal 
crossover between weak- and strong-correlation regimes. In this connection the ab- 
sence of Brinkman-Rice transition, which arises when one resorts to an additional 
mean-field-type approximation (so-called Gutzwiller approximation) in estimating 
expectation values, is checked again based on the VMC data in ID, 2D and 3D 
with better statistics in Appendix A. Furthermore, in Appendix B we confirm the 
absence of long-range orders within GWF. In §4 to improve the normal state on 
GWF we introduce a binding factor between j- and J,-spin electrons by taking ac- 
count of perturbation processes in the strong coupling limit. This wave function not 
only remedies the shortcomings of GWF but undergoes a substantial metal-insulator 
transition. As a supplement to §4 details of the metal-insulator transition in this im- 
proved wave function are described in Appendix D. As an alternative improvement 
on GWF an distance-dependent intersite correlation (we call 'Jastrow') factor are 
studied in Appendix C. Section 5 is assigned to summary. 

§2. Model and Method 

2.1. Attractive Hubbard model 

In this paper we study a couple of basic wave functions as a normal state for the 
attractive Hubbard model (U < 0): 

H = H t + Hu = -t ]T (c\ a c ha + c) a a, a ) + Uj2n jmi , (2-1) 

where n Jcr = c^ a Cj a and other notations are standard. We use t as the unit of energy. 
In this paper only bipartite lattices are taken up, and the sum of transfer in eq. (2-1) 
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is restricted to the nearest-neighbor pairs. We concentrate on the 2D square lattice 
except for Appendix A, where the dependence on lattice dimension becomes an 
important element. Due to the electron-hole symmetry we restrict electron density 
n (= N c /N) to n < 1 with N c and N being the numbers of electrons and sites, 
respectively. Here, we exclude an external field, thus iVf = iVj_ = N e /2 (N a : number 
of cr-spin electrons). A chemical potential term —QYl<ja n j<T ma Y De added to adjust 
electron density, if necessary. For later use we summarize a couple of basic facts 
about this model in the following. 

For bipartite lattices, AHM is related to RHM by a celebrated canonical trans- 
formation, n )' 12 ) 

cj1=b jh c n = e iG ^b) v (2-2) 

where G is the antiferro reciprocal lattice vector, e.g. G = (vr,7r) for the 2D square 
lattice, and rj is the position of site j. The transformed Hamiltonian is written as, 

K = -*E + 6 5 A*) - U £ h n h n + LWn T - h ]T fe + I) , (2-3) 

where = bj a bj a , S*| = (n^f — hj^)/2 and n CT = N a /N. Henceforth, representation 
with a tilde denotes the transformation by eq. (2-2). Thus, the sign of the inter- 
action term is reversed, and the chemical potential £ (electron density n) in AHM 
is related to the effective magnetic field as h = 2( (magnetization as m = 1 — n) 
in z direction in RHM. Moreover, unless the original AHM possesses spin polariza- 
tion, electron density in the transformed RHM is always at half filling due to the 
relation, n,j + = 1 + rijj — jijj,. Meanwhile, the order parameters of antiferro 

CDW and onsite singlet pairing defined as, Ocdw = j? ^2j(~ 1) J ^( n jT + n jl ~ ^) 

and Osc = jj J2j( c ]-\ c ]i) or jjJ2j( c jl c jl)7 are transformed into the forms of z and 

x,y components of SDW order parameter, Og D w = ~N l) J ( n jT ~ n jl) an< ^ 

^sbw = ^Sj( c ]| c j|) or ~k ( c ]| c it)' respectively. 13 -* Through this transforma- 
tion one can deduce the properties of AHM from the knowledge of RHM, and vice 
versa. In particular, at half filling (( = 0) the relations become simple owing to the 
disappearance of the effective field. For example, it is widely accepted that an antifer- 
romagnetic long-range order with a common magnitude of Ogrjw amon g a = x iVi z 
arises in the half-filled RHM on 2D square lattice for arbitrary U (> 0). This reads 
as AHM at half filling simultaneously possesses a CDW and a singlet pairing order of 
the same magnitude. 13 ) Note that it is restricted to exact treatments that the above 
mapping holds unconditionally; when some approximation is applied, the validity of 
mapping has to be verified against individual treatments. 

Next, let us consider the strong coupling case of eq. (2-1). As the t-J model is 
derived from RHM, an effective Hamiltonian for large- \U\/t regime is obtained by 
picking out the lowest-order terms in a perturbation expansion in t/\U\ as, 14 )' 15 ) 

fteff = TfJT {-PiPj + PiPj) > ( 24 ) 
11 («> 
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where pi = c^Cii and pi [= {n,^ + n^)/2] is the corresponding number operator. 
They operate in the subspace with no singly-occupied site. Thus, the system is 
reduced to hard-core bosons with nearest-neighbor repulsive interaction. Notice that 
in the strong-correlation region where eq.(2-4) is valid, every energy within 'p'-band 
is scaled by t 2 /\U\, which is much smaller than the energy of order \U\ needed to 
pull a pair apart: the Hubbard gap. 

2.2. Variational Monte Carlo method 

We use a conventional variational Monte Carlo (VMC) method. 16 ) However, 
the accuracy of calculations in this paper is notably improved by virtue of the rapid 
progress of computer technology in the past decade. In this work we chiefly consider 
the 2D square lattice with N = L x L sites, and take up the 3D simple cubic lattice 
with N = Lx Lx L only in Appendix A to see space-dimension dependence of GWF. 
To discuss the thermodynamic limit we carefully check the system-size dependence 
with systems up to N = 2N C = 28 x 28 = 784 and N = N c = 24 x 24 = 576 for 2D 
and N = N c = 10 x 10 x 10 = 1000 for 3D. 

The periodic(P)-antiperiodic(A) boundary condition is imposed on all the sys- 
tems we use for the square lattice. Thereby, systems with arbitrary even L for n = 1 
can be treated satisfying the closed shell structure, which enables us to check system- 
size dependence. On the other hand, for n = 0.5 available systems are limited to 
L = 8, 12, 20, 28 etc. under this condition, but we may regard the system of L = 28 as 
sufficiently large as will be mentioned later. For incommensurate filling it is practi- 
cally impossible to directly check system-size dependence, because available densities 
are different size by size. For the 3D simple cubic lattice in n = 1 the P-P-A and 
the P-A-A boundary conditions are used for lattices with L = 41 and L = 41 + 2 (/: 
integer), respectively, to satisfy the closed shell condition. To reduce the statistical 
error a great number of samples are used: 2.5 x 10 5 to 10 according to N e . To keep 
the mutual independence in electron configurations among samples, three kinds of 
trials (single electron hopping, exchange and pair hopping) are adopted; besides, the 
sampling intervals are adjusted up to 100 Monte-Carlo steps so that each acceptance 
ratio may reach 50%, except for some limiting cases. 

§3. Gutzwiller wave function 

In this section as a continuation of the previous researches for RHM 16 ) we study 
the properties of the Gutzwiller wave function for AHM. In §3.1 we briefly summarize 
fundamentals of GWF for later use. In §3.2 and §3.3 we argue the VMC results for 
energy and correlation functions, respectively. 

3.1. General aspects 

As a trial wave function for normal phases of AHM, the celebrated GWF, 17 ) 
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is again a useful starting point due to its simplicity. In eq. (3-1), g or a(= Ing) is the 
sole variational parameter, which controls the ratio of doubly occupied sites d, and 
<Pf is the Fermi sea. For AHM the range of g will be 1 < g < oo to enhance d, whereas 
for RHM < g < 1 to reduce d. The two limiting cases, g = and g = oo, indicate 
insulating states without doubly-occupied sites and without singly-occupied sites, 
respectively. The former state is used for the t-J model, 18 ) and becomes the exact 
ground state particularly for the ID super symmetric model with 1/r 2 coupling. 19 ) 
For g = 1 GWF is reduced to the noninteracting case, <^f- 

So far, GWF has been extensively studied for RHM to treat itinerant mag- 
netism, Mott transition, 3 He, superconductivity etc. Because of the difficulty in 
calculating the expectation values by GWF, they had been estimated by a mean- 
field-type approximation [Gutzwiller approximation (GA)] for a long time since the 
first introduction by Gutzwiller himself. 20 )' 21 ) Such an additional approximation, 
however, breaks the variation principle and sometimes leads to qualitatively wrong 
results. A typical erroneous case is the Brinkman-Rice metal-insulator transition, 22 ) 
which never takes place actually in finite dimensions. We will reconsider this point 
in detail in Appendix A. Afterward, accurate calculations have been carried out by 
a VMC method 16 ^ and diagrammatic calculations for D = 1 and oo. 23 )' 24 ) As for 
AHM, the behavior of GWF has not yet been addressed sufficiently, although GA 
has been already applied. 25 ) 

To take advantage of the fruits of RHM, let us look into the relationship between 
GWF's in the two models through the canonical transformation eq. (2-2), first. Since 
the relations of operators in k space are = and Cki = b^_ k+G ,, the vacuum 

state for b^ operators is defined as |0)| = life c ijO)l and inversely |0) = life &fcj_|0)l- 
Then, the transformed Fermi sea is written as, 

n b U n h + GiU b i>ifi)- ( 3 - 2 ) 

fc<fcp fc<fcp k 1 

Except for half filling, <Pp has a finite magnetization according to the effective 
field h in eq. (2-3). In particular, if the nesting condition is completely satisfied, 
e.g. for the square lattice and the simple cubic lattice at half filling, <Pf becomes 
equivalent to <?f- On the other hand, the Gutzwiller factor is transformed into 
Vq = exp (aJVj) exp (^—aJ2j n j^jlj j due to (J2j n j|)^F = N^<Pp. Consequently, by 
neglecting the constant factor, the transformed GWF is written as 



*g = n 1 - ( x - 1) 



<P F . (3-3) 



By comparing I^g with the original GWF eq. (3-1), the Gutzwiller parameter g 
is replaced by l/g{= 7), in addition to a reversal of the Fermi sea eq. (3-2). In 
particular, if the system has the electron-hole symmetry, the relation, 

E ti9 ) = Et(7 ), ^2>_i_*M (34) 

is derived from the equivalence between $g(<7) and <Pg(7) and eq. (2-3). Here, Et{g) 
and Ejj(g) are the variational expectation values of TL t and Jiy with respect to 



6 



Hisatoshi Yokoyama 



\pG(g), respectively. Thus, the energetics for RHM is simply mapped to the one for 
AHM in this special case. In Figs. 1(a) and (b), E t and Ejj by GWF for the square 
lattice are shown as a function of log 7. The symmetrical feature with respect to 
7 = 1/(? = 1 given by eq. (3-4) is observed for half filling in both components. As n 
decreases from 1, this symmetry is deservedly broken. Henceforth, we will often use 
7 instead of g, because 7 in AHM plays a substantially equivalent role to g in RHM. 

Fig. 1. Expectation values of (a) Tit and (b) Hu with respect to GWF eq. (3-1) are shown for 
several electron densities. The data for 7 < 1 (7 > 1) correspond to the attractive (repulsive) 
case. Samples as many as 5 x 10 5 -5 x 10 6 are used for each data point in VMC calculations. 
Statistical errors are negligible in this scale. 

3.2. Results of energy 

Now, we reconsider the energy of GWF by VMC calculations. In Fig. 2 expec- 
tation values of kinetic and potential energies are plotted versus 7 in the attractive 
regime. Since onsite pairs are more favored as 7 decreases from 1, the ratio of double 
occupation d (= Ejj/U) increases, and the system loses kinetic energy. In the limit 
of 7 — > 0, all the electrons form onsite pairs for arbitrary electron density: d = n/2. 
Thus, Et vanishes for any value of n. This feature is in sharp contrast to the repul- 
sive case, in which Et has finite transfer even at g = except for n = 1. Since the 
behavior of Et/t and Ejj/U in this limit is essential for the Brinkman-Rice metal 
insulator transition (BRT), we will investigate their behavior in detail in Appendix 
A. 

Fig. 2. Expectation values of (a) Ht and (b) Hu by GWF for all the possible electron densities 
for L — 10 in the attractive regime. 5 x 10 5 -10 7 samples are used for each value of 7 in VMC 
calculations. Statistical errors are by far smaller than dots. 

Total energy E, the sum of E t and Ejj, has a minimum at < 7 < 1 for any finite 
value of negative U, because BRT never arises in finite dimensions as discussed in 
Appendix A. As an example Fig. 3 shows E/t for two values of U/t at half filling; the 
values for L = 00 (cross) are estimated by some extrapolations. 26 ^ With increasing 
\U\/t the optimal value of 7 naturally decreases. Here, note the difference between 
Figs. 3(a) and (b) that the optimal value of 7 hardly depends on L for U/t = —8, 
while it remarkably decreases with increasing L for U/t = —16. By checking the 
data for various values of U/t, it is found that the former (latter) feature is common 
to the regime of small (large) \U\/t, and a rather precipitous crossover takes place 
at U = U co lit. 

Such a crossover is recognizable also in Fig. 4(a), where the optimized values of 
7, E t /t and Ejj/U are shown for a couple of system sizes. With increasing \U\/t, 
7 at first shows a rapid decrease almost without system-size dependence. When U 
reaches U co , however, the decrease of 7 becomes slow and asymptotic with respect to 
7 = 0, and simultaneously system-size dependence becomes conspicuous. With this 
change of 7, \E t \/t decreases with weak system-size dependence originating in the 
noninteracting band (not by 7), whereas large dependence shows up for \U\ ^ \U co \ 
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Fig. 3. Total energy of GWF at (a) U/t = —8 and (b) —16 for several system sizes at half filling. 
Symbols are common in (a) and (b) . Minimal points are indicated by solid circles connected by 
dashed lines. The crosses on these lines indicate the values in the thermodynamic limit, which 
are extrapolated from the data of every even L (6 < L < 20) and also L — 24 for U/t = — 16 
by second-order polynomial fit of 1/L 2 . The value of GA is also shown in (a); as for (b), the 
total energy of GA is U/2 = -8t because \U\/t > \U B R.\/t = 12.969 • • •. Incidentally, although 
E/t of GWF is always lower than the value of GA in Fig. 3, this is not always true, because GA 
violates the variation principle. 16 - ) 



due to the behavior of 7. As for the potential part, double occupancy d increases 
and almost saturates for \U\ ^ \U co \, and its system-size dependence is reversed at 
U~U co . 

Fig. 4. Optimized values of 7 and the energy components Et/t and Eu/U as a function of \U\/t at 
(a) half filling and (b) a quarter-filling. For clarity only a couple of system sizes are taken up. 
The results of GA are also plotted, in which BRT takes place at Ubb.- 

Plotted in Fig. 5 is the total energy E/t measured from that of the completely 
localized or paired state: Eoo = Un/2. The magnitude of this quantity rapidly 
reduced with increasing \U\/t for small \U\; around U = U co , however, the behavior 
is switched to —t/\U\. Thus, GWF represents two kinds of normal states, according 
as \U\ ;€ \U co \ or \U\ ^ \U co \. In the former region the state obviously belongs to 
a Fermi liquid, which continues to the noninteracting case, whereas in the latter 
most of the electrons form onsite pairs as in Fig. 4(a), and the state looks consistent 
with the picture of hard-core bosons represented by eq. (2-4), in which the kinetic 
contribution does not vanish but is proportional to t 2 /U. As will be seen in the next 
subsection, however, this state is still metallic — an almost localized Fermi liquid — , 
because a clear Fermi surface is defined and there is no excitation gap. This is in 
contrast with the GA [solid line in Fig. 4(a)], which brings about a metal-insulator 
transition at U = Ubr- In the state for \U\ > \Ubr\ all the electrons stand still and 
Ef vanishes, although it is quantitatively good for \U\ < \Ubr\- This unfavorable 
feature is the same with the insulating state obtained by DMFT. 7 )' 8 ) 

Fig. 5. Total energy of GWF at half filling as a function of \U\/t measured from E^ — nU/2, 
namely the energy of the completely paired state. Shown is the value for L — 00 extrapolated 
from the data of L = 6-20 and 24 (for large values of \U\/t) as in Fig. 3. The thin dashed line 
with an arrow shows the — t/\U\ curve whose coefficient is chosen so as to fit the data of L — 00. 
The inset shows the magnification of small- E region; the data for some finite L's are included. 
Incidentally, this quantity is equivalent to the simple total energy of RHM at half filling [See 
eq. (3-4)]. 

Now, we turn to less-than-half filling. Plotted in Fig. 4(b) is the same quantities 
with Fig. 4(a) but for a quarter filling. As compared with half filling, the qualitative 
feature scarcely changes including system-size dependence. Therefore, fixing the 
system size, we focus on the n dependence. In Fig. 6 we show the optimal values of 
the energy components versus \U\/t for every available value of n for L = 10. Both 
E t and Ejj change their behavior around U = U co from precipitous increase with 
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increasing \U\/t (U ^ U co ) to almost saturated one (U ^ U co ). As has repeatedly 
insisted, in the former region the kinetic term is dominant (Fermi liquid), while in 
the latter region the potential term is superior (almost localized). An important 
point here is that although the absolute values of energies are directly related to the 
electron density, the values of |f7 co | once slightly increases with decreasing n, has a 
broad maximum near n = 0.5 and shifts to a somewhat smaller value in the vicinity 
of n = 0. On the whole this behavior agrees with that of |C/br| in GA (Fig. 19), 
although \U co \ is smaller than |C/br| to some extent. To summarize, in contrast to 
RHM, an almost localized state exists also for less-than-half filling in a similar range 
of U/t to half filling. 

Fig. 6. Energy components (a) E t /t and (b) Eu /U (= d) at the optimal 7 as a function of \U\/t 
for all the electron densities available for L — 10. The results of <Pq (L = 10, n = 1.0) treated in 
§4 are also shown for comparison by shaded dashed lines. Symbols are common to both figures. 
Energy minimization is carried out based on the data points with 10 6 -5 x 10 6 samples. 

3.3. Results of correlation functions 

In this subsection we study GWF through momentum distribution, n(k) = 
I J2a( c ka c kcr): an d a couple of correlation functions defined as, 

S z (q) = ^Y1 e ~ iqn ( S j S j+^> ( s P in z component) (3-5) 

jl 

N(q) = ^ E e " iqn ( n i n J+e) ~ n2 > ( char S e density) (3-6) 

jl 

ji 

where rij = nji + n,| and Aj = cj^c^. Some properties of these quantities have 

been already studied with ^ G for the ID and 2D RHM. 16 )' 27 ) Although they partly 
exhibit unphysical behavior, 28 ) they are suggestive and useful to compare AHM with 
RHM. 

To apply the previous results of RHM to AHM, we summarize the relationship 
among correlation functions through the canonical transformation eq. (2-2). Since 
N t = (n jnj+e ) - {nj)(n j+t ) is transformed to Sf = 4((S*S] +t ) - (S*)(S* +t )), N(q) in 
AHM corresponds to S z (q) in RHM, as anticipated from the relation between Oqdw 
and Ogrjw Similarly, P(q) in AHM is related to the staggered spin correlation in 
RHM as, P(q) = ^(q + G)/2, because A]A j+t , is transformed to (-l) j Sj~S~ +e , 
where 

S± (1) = ^E-- iqn (StS- +e + S-Sl +e ). (3-8) 
jl 

Using these relations and assuming S{q) = S z {q) = ^(q), we have N(q) = 2P(q — 
G). In particular, the case for q = G [N(G) = 2P(0)] tells that the correlation of 
antiferro CDW is always proportional to the correlation of onsite singlet pairing, if 
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the wave function, which may be exact or approximate, remains isotropic in spin 
space after the transformation. Since this condition is satisfied by GWF at half 
filling due to eqs. (3-2) and (3-3), the relation [S(q)] g = [N(q)] 7 = 2[P(q - G)] 7 
holds, where [• • -] g indicates the expectation value with respect to W r Q(g). However, 
this is not the case for n ^ 1, where does not meet the condition for spin isotropy. 
As for momentum distribution, particularly at half filling one can confirm the relation 
[n(k)] g = [n(/c)] 7 by using the inversion and electron-hole symmetries in addition to 
the transformation eq. (2-2). 

To begin with, we look at the half-filled case for U < 0, first. In Figs. 7(a) we 
show n(k) for several attractive values of g, which are equivalent to the repulsive case 
with 1/g or with \U\/t. Note that the discontinuity at k = kp is clear for \U\ > \U co \, 
not to mention \U\ < \U co \. In Fig. 17 we show the renormalization factor Z (quasi 
particle weight in the Fermi liquid theory) estimated from the discontinuity of n(k) 
at k = k?: Z = n(/cF — 0) — n(/cF + 0); Z of GWF decreases similarly to GA for small 
values of \U\/t, but abruptly changes the behavior around U = U co ~ —lit to an 
asymptotic one. Thus, we confirm GWF remains metallic also in the strong-coupling 
side of Uco, even though the effective mass of electrons (1/Z) is considerably large. 

Next, we touch on correlation functions. According to the previous studies for 
U > 0, N(G) and P(0) should be enhanced and S(q) is suppressed by attractive 
correlation; this aspect is actually seen in Fig. 7(b), where N(G) becomes over- 
whelmingly large. Therefore, we check the possibility of the divergence of N(G) as 
well as the realization of a long-range order in Appendix B. Here, we focus on the 
gaps in excitation spectra, which we learn from the behavior of N(q) and S(q) near 
the T point (\q\ — ► 0). In this limit of |g|, both N(q) [inset of Fig. 7(b)] and S(q) 
[N(q) in Fig. 1(a) of ref.27)] behave as a linear function of \q\ irrespective of the 
values of U/t. Therefore, GWF is always gapless both in spin and charge degrees of 
freedom, supporting GWF's metallic properties in both sides of U co . This aspect is 
preserved in less-than-half filling as observed in Figs. 10(a) and (d). 

Fig. 7. (a) Momentum distribution function and (b) charge density structure factor of GWF at 
half filling for several attractive values of g; the corresponding values of U/t [common to (a)] are 
shown in (b). The values for \U\ < \U co \ {\U\ > \U co \) are shown by solid (open) symbols. Inset 
in (b) shows the magnification of N(q) on the T-X line. In (a) each k point shifts by tt/L in y 
direction due to the antiperiodic boundary condition; a mark with a prime (e.g. F') indicates a 
shifted point. The arrow shows the position of M' point. 

Now, we turn to the less-than-half filling, where the relation to RHM is not so 
simple as at half filling, and N(G) and 2P(0) are no longer equivalent. Since the 
correlations of charge density and onsite singlet pairing are dominant in AHM, we 
take up N(q) and P(q), first. To grasp the n dependence we fix the parameters 
at U/t = —8 and L = 10 (Fig. 8). With decreasing n the sharp peak of N(q) at 
q = G is rapidly suppressed, because the nesting by G is lost. Instead, weak singular 
behavior appears at incommensurate wave numbers [q = 2kp or 2(G — kp)] of the 
noninteracting case for n ^ 0.72. 29 > When n further decreases, N(q) tends to be small 
and flat, namely charge correlation becomes weak and restricted to a considerably 
short range. We show \U\/t dependence of N(G) for L = 10 in Fig. 9(a). At half 
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filling, N(G) shows precipitous increase around U = U co , but slowly saturates for 
\U\ ^ \U co \. When n decreases and G (M point) is no longer a special g-point, the 
conspicuous change around U ~ U co vanishes. 

In contrast to N(q), P(q) preserves the peak at q = (0,0) even for the lowest 
electron density n = 0.04 as shown in Fig.8(b), although its magnitude gradually 
diminishes with n. Thus, the singlet pairing correlation remains dominant for all 
n in GWF, although it does not develop into a long-range order (see Appendix B). 
As for \U\/t dependence [Fig. 9(b)], P(0) abruptly increases around U = U co for 
every electron density. Note that for \U\ > \U co \ the absolute values of P(0) scarcely 
changes with n for n ^ 0.72, and has rather larger values for n = 0.96 and 0.88 than 
for half filling. Concerning S(q), which is considerably suppressed and flattened by U 
already at half filling, the tendency toward shortened correlation is further promoted 
by decreasing n [inset of Fig. 8(b)]. 

Fig. 8. Structure factors of (a) charge density and (b) onsite singlet pairing and spin (inset), for 
all the electron densities available in L = 10. Symbols are common to all the panels. In each 
panel the value for U/t = (L = oo) at half filling is shown by a dashed line. We take a path 
r^X^M->Tina Brillouin zone. 5 x 10 5 (n = 1.0)-5 x 10 6 (n = 0.04) samples are used. 

Fig. 9. Structure factor of (a) charge density at q — G and (b) onsite single pairing at q = as a 
function \U\/t for all the electron densities available in L — 10. 

Finally, we consider the difference between AHM and RHM for n < 1. Shown 
in Figs. 10(c) and (d) are n(k) and S(q) of GWF at a quarter filling for several 
attractive values of g, respectively. One finds similar features to those at half filling, 
for instance, [1] n{k) is an increasing function of and [2] has clear discontinuities 
at k = kp for any finite value of g. [3] Z changes its behavior at U ~ U co (See Fig. 17). 
[4] S(q) is remarkably suppressed with increasing correlation strength \U\/t. 

As for [1], the increase of n(k) especially inside the Fermi surface in Fig. 10(c) is 
not only common to but more remarkable than that at half filling [Fig. 7(a)]. This 
tendency is contrary to that in less-than-half filling of RHM [cf. Fig. 2(b) of ref. 27)]. 
28 ) A similar situation, that is, the feature is analogous between n < 1 and n = 1, is 
also observed in N(q) [Fig. 10(a)] and S(q) for AHM. As mentioned, N(q) comes to 
have higher peaks (though not so dominant as for n = 1 [Fig. 7(b)]) at k = 2kp and 
2(G — kp) with increasing g. Meanwhile, S(q) in less-than-half filling for RHM does 
not have peaks at those wave numbers but only cusps [See Fig. 2(a) of ref. 27)]. As 
for S(q) — the point [4] [Fig. 10(d)] — , the maximum for g — 1 = 50 (U/t ~ —15) is no 
more than 1.6% of the maximum value 0.5 for g — 1 = (U/t = 0); the corresponding 
value at half filling [g- 1 = 50 (U/t ~ -15), L = 20] is 0.85%. In both densities, the 
reduction of S(q) is considerable. On the other hand, the corresponding quantity 
[N(q) for U/t = 16, L = 16] at a quarter filling of RHM preserves its amplitude as 
large as more than 60% of the free case [Fig. 2(a) of ref. 27)]. 

Accordingly, the state for n < 1 in AHM has analogous properties of correlation 
functions to its counterpart in RHM at half filling, but not for n < 1 in RHM at 
least within GWF. 
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Fig. 10. Structure factors of (a) charge density, (b) onsite singlet pairing and (d) spin, and (c) 
momentum distribution by GWF at a quarter filling for several values of attractive g (7), which 
correspond to the values of U/t shown in (b) and (c). Arrows on the noninteracting case (g = 1) 
indicate the points of 2fcp or 2(G — far). The corresponding singular behavior in N(q) is indicated 
by the letters A, B and C in (a). Path in q space is the same with Figs. 8 and 7(a). 2.5 x 10 5 
samples are used for each parameter set. 

§4. Intersite correlation factor 

As discussed in our previous papers 16 )' 27 ) GWF is too simplistic to reproduce 
properly some aspects of momentum distribution and correlation functions. 28 ^ This 
is primarily because intersite correlation factors are neglected in GWF, as we showed 
for RHM. In this section and Appendix C we improve the normal-state wave function 
on GWF by taking account of these factors. In Appendix C we take up a distance- 
dependent Jastrow-type correlation factor, which was successful for less-than-half 
filling of RHM. 27 ) On the other hand, in this section a binding effect between up 
and down spins in adjacent sites is introduced by allowing for the virtual processes 
in the perturbation expansion in the strong coupling limit. In §4.1 we introduce the 
binding wave function $q. In §4.2 and §4.3 the results at half filling and in less- 
than-half filling are presented, respectively. Section 4.4 is assigned to discussions 
and comparisons with other results. A further discussion on the metal-insulator 
transition in and beyond "Pq is given in Appendix D. 

4.1. A binding factor between up- and down-spin sites 

Consider the strong coupling limit (t/\U\ — > 0). All the electrons form onsite 
pairs in the ground state [Fig. 11(a)]. A small transfer term causes an on-site pair 
broken, yielding thereby virtually excited |- and J,-spin pairs next to each other 
[Fig. 11(b)]; they corresponds to the virtual states in the lowest-order perturbation 
processes. Such a configuration should have a larger weight than ones for higher- 
order processes, which include isolated singly-occupied sites as in Fig. 11(c). In fact, 
improvement in this line has been already embodied for the half-filled RHM, 30 )' 27 ) 
where the binding between doubly-occupied (d-) and empty (e-) sites is vital to the 
virtual processes. The wave function thus introduced (<^q) 27 ) was extremely good, 
especially for the ID half-filled band, not only in energy but for various correlation 
functions. 

Fig. 11. Schematic examples to assign weights of the correlation factor in \Pq on 2D square lattice. 
Each figure represents a local electron configuration, where an open circle indicates an empty 
site, a solid circle a doubly-occupied one and an upward (a downward) arrow a singly occupied 
one with an up (a down) spin, (a) A configuration realized in the ground state at U/t = — 00. 
Here, we assume the weight of this configuration as 1. (b) A configuration in which one electron 
pair [marked with a cross in (a)] is broken by a one-step hopping process (a virtual state in the 
second-order perturbation processes). The weight 7(= 1/g) is assigned due to the Gutzwiller 
factor, (c) A configuration in which one electron pair is broken and two hopping processes are 
needed to reach from (a) (a virtual state of the fourth-order perturbation processes). The weight 
is suppressed by the factor (1 — fi) 2 in addition to 7. 
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For the attractive case a suitable wave function is expected by applying the 
canonical transformation eq. (2-2) to $q, and consequently we have, 

^0 = I1[ 1 -^]^G(7), (4-1) 

j 

^-^na-^+^IK 1 -^)' (4-2) 

T T 

where sj = rij a (l — nj- a ) and r runs over all the nearest neighbors. In \Pq the 
second parameter fj, (0 < ji < 1) controls the attractive correlation between a pair of 
up- and down-spin electrons in next-nearest-neighbor sites both of which are singly 
occupied as explained in Fig. 11. When /j, = 0, $q is reduced to GWF. As \i increases, 
singly-occupied sites without an antiparallel-spin site in nearest neighbors become 
disadvantageous. The state becomes insulating at last for fx = 1, because then an up- 
spin and a down-spin electron are completely bound within nearest-neighbor sites, 
unless they are on the same site. A merit of among similar functions 31 ) is that the 
on-site correlation and the nearest-neighbor binding can be treated independently 
by g and fi, respectively. 

In particular, at half filling the relation eq. (3-4) holds also between $q(<?, /i) and 
$q(7, h), because the correlation factor for binding in eq. (4-1) is invariant under the 
canonical transformation and GWF satisfies eq. (3-3). Thus, the results obtained in 
RHM 27 ) can be directly applied to AHM at half filling. For example, is lower 
in energy than an antiferromagnetic state (i^sdw) 32 ^ m ID, but higher in 2D for 
RHM. 27 ) This reads for AHM as <Pq is stabler than "Pcdw 33) in ID, but not in 2D. 

4.2. Results for half-filled band 

First of all, let us reconsider the case of half filling. In Fig. 12 expectation values 
of the energy components are shown for various values of the parameters. As for 
the kinetic part, Et once slightly decreases with n and have a minimum at a small 
value of \i for 7 < 1 [inset of Fig. 12(a)]. This decrease of E t causes the gain in 
total energy for small \U\. For /x ^ 0.2, E t comes to increase monotonically, but 
remains finite even for \i = 1 due to the second-order perturbation processes. On 
the other hand, double occupancy Ejj/U decreases with increasing [i for 7 0.5 
due to the fluctuation effect of /x, whereas d increases with fj, for smaller values of 7 
[inset of Fig. 12(b)], implying that isolated unpaired electrons are confined in local 
singlet pairs. Anyway, $q has lower energy than GWF for any negative value of U, 
as we showed in the previous study for RHM. The energy decrement is remarkable 
especially in the intermediate correlation regime [See Fig. 22]. In Fig. 13 total energy 
for various parameter values are depicted as an example for large- \U\ region. 

Fig. 12. Expectation values ol (a) kinetic and (b) potential energies of •Z'q as a function of 7 = 1/g 
at half filling. Insets show /1 dependence for respective quantities. Symbols are common in (a) 
and (b). Although the system with L = 6 is plotted for abundant data points, it is confirmed 
that the system-size dependence is quantitatively small. We use 10 6 samples for each data point. 
Statistical errors are much smaller than the symbol size. 
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Fig. 13. Example of energy expectation values of <Pq for a large value of \U\/t as a function of 
binding correlation strength fi for several values of onsite attractive correlation g. The arrow 
indicates the minimum -8.2118(23), which is given by [i = 0.9 and 1 — g = —7.0 and much lower 
than the value of GWF (fi = 0) -8.0463(12) but comparable to the insulating (/i = 1) value 
-8.2115(23). 

Next, we look at the optimized variational parameters (Fig. 14). For small 
\U\( ^ W) the optimized value of [i greatly increases with increasing \U\/t, and its 
system-size dependence is small. For \U\ ^ W, however, \x comes to change slowly 
near 1, and becomes highly dependent on L and n. On the other hand, the optimal 
value of 7 for I^q is almost the same with (slightly smaller than) that of GWF for 
\U\ ^ W, undergoes a cusp- like change at U ~ W and becomes considerably larger 
for \U\ ^ W. The dependence on L and n is weak for any value of U. These 
facts indicate that there is some crossover at \U\ = \Uq\ ~ W and the state in the 
strong-coupling side is qualitatively different from GWF, and nearly insulating. 

Fig. 14. Optimized values of the parameters in versus \U\/t for some system sizes and densities. 
Conspicuous errors in /i for \U\ il \Uq\ stem from the flat bottom of E/t near /i = 1 (Sec §4.4). 
For comparison, the optimal values of 7 for GWF are also shown. 

To study this crossover closely we consider energy components E t /t and d, which 
are shown in Figs. 15(a) and (b), respectively. With increasing \U\/t, \Et\ abruptly 
decreases until \U\ reaches \Uq\, around which \Et\ changes its behavior, and comes 
to decrease slowly for \U\ > \Uq\. Note that as L increases, the change of behavior 
at U ~ Uq becomes more abrupt and cusp-like. A similar abrupt change at U ~ Uq 
can be seen in Ey. Thus, an abrupt crossover or a phase transition 34 ) takes place 
around U = Uq, namely I^q represents two distinct normal (or incoherent) states. 
The critical value Uq is broadly estimated as —8.8 at half filling by finding the 
intersection of a smoothly extrapolated curve from \U\ < \Uq\ and a line of fi = 1 
with respect to d [Fig. 15(b)], which we adopt on account of its sharper change 
there. The value thus obtained does not largely depend on L. As will be explained 
below with correlation functions, in the region of \U\ < \Uq\ the state is a Fermi 
liquid, which is smoothly connected to the noninteracting case, whereas the state for 
\U\ > \Uq\ is very akin to the insulating state (fi = 1) indicated by dashed lines in 
Figs. 15. In this nearly insulating state almost all the electrons are paired either at 
the same site or within the nearest-neighbor sites. As compared with the optimal 
GWF (Fig.6), double occupancy in 'J'q has a larger value around U = Uq but is only 
a little smaller for \U\ ^ \U co \(~ lit) due to the fluctuation effect of fi. Meanwhile 
!Pq is considerably advantageous in kinetic energy in the strong coupling side, thanks 
to the second-order perturbation term. 

Now, we look at this transition in momentum distribution function shown in 
Fig. 16 for various values of U. In the region of |£7| < \Uq\ discontinuities at the two 
kp points are easily recognized, indicating a degenerate metallic state. On the other 
hand, for \U\ ^ \Uq\ the discontinuity at kp vanishes or is at least not obvious, so 
that the state is insulating. Furthermore, near this boundary there is a qualitative 
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Fig. 15. Expectation values of (a) kinetic and (b) potential energies by the optimized parameter 
sets as a function of \U\/t at half filling. Four system sizes are shown. Dashed lines denote the 
optimized insulating values (/i = 1). Insets show the same quantities for four values of electron 
density with a fixed system size L = 10. We have used 5 x 10 5 -2 x 10 6 samples for energy 
minimization (each parameter set). 

distinction in system-size dependence between the two region (not shown), namely 
the discontinuity at kp increases with increasing L for \U\ < \Uq\ (e.g. U = —8t), 
whereas it rapidly decreases and tends to vanish for \U\ > \Uq\ (e.g. U = —9t). In 
Fig. 17 renormalization factor Z is depicted (solid circle) as a function of \U\/t. Z 
comes to decrease remarkably as U approaches Uq and vanishes there; the effective 
mass of quasi particles diverges and the state becomes insulating. Thus, the existence 
of a metal-insulator transition is ascertained (See Appendix D for details). 

Fig. 16. Momentum distribution function by at half filling for various values of U/t. Closed 
(open) symbols are used for the Fermi liquid (insulating) regime; U/t = —9 is near the critical 
point. The path in q space is the same with Fig. 7(a). We use 10 6 samples for each optimized 
parameter set. 

Fig. 17. Renormalization factor (quasi particle weight) in Fermi liquid states versus \U\/t estimated 
from the discontinuity of n(k) at k = fcp on the T-X line for some electron densities. Results of 
GWF and GA in 2D and of DMFT for hypercubic lattice 35 ' are also included. 

Let us move on to N(q), 36 ^ which has a characteristic peak of half filling at the 
M point as seen in Fig.18. In the inset N(G) is plotted versus U/t; the increase of 
N(G) with \U\/t is abrupt right below \U\ = \Uq\, but becomes slow for \U\ > \Uq\. 
Thus, the effect of transition is clear also in N(q). At these sizes the magnitude of 
N(G) for large \U\ is about an order of magnitude smaller than that of $j in the 
CDW phase (Fig. ??), is not proportional to system size even though increasing with 
L, and is the same degree with that of GWF (Fig. 9). In addition, Ocdw is always 
zero irrespective of the value of L within the statistical error. It follows that only 
short-range correlations grow even for \U\ ^ \Uq\ in 

Finally, we look at excitation gaps. As for the charge degree of freedom, N(q) ~ 
v\q\ for \q\ — > as seen in Fig. 18, and v only weakly depends on the value of \U\/t. 
Thus, there is no charge gap in both phases. On the other hand, there is a distinction 
in the small-|g| behavior of S(q) between the metallic and the insulating phases, as 
observed in the inset of Fig. 20(b). For \U\ < \Uq\ the leading term of S(q) seems 
linear (or a smaller exponent), while for \U\ > \Uq\ it seems S(q) tx \q\ a with a > 1. 
Thus, we conclude that a spin gap does not exist for the metallic case, but does for 
the insulating case, in which finite energy is necessary to break a local singlet pair. 

Fig. 18. Charge density structure factor by tf'Q at half filling for various values of U/t. Inset shows 
N(q) at q — G (M point) as a function of coupling strength. Four system sizes are shown. 
10 6 -2 x 10 6 samples are used. 
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4.3. Results for less-than-half filling 

In the repulsive case $q becomes inadequate as soon as the electron density 
leaves half filling, because the symmetry between d-sites and e-sites breaks. In the 
attractive case, however, !Pq remains good even for less-than-half filling as will be 
shown, because f-site and j-site preserve the symmetry even for n < 1. 

First, we touch on the behavior of energy as a function of [i (not shown). For 
n ^ 0.5 the feature of both E t /t and Ejj /U is qualitatively the same with at half 
filling (Insets of Fig. 12), whereas for densities as low as n = 0.24 both E t /t and 
Ejj/U are always monotonically increasing functions of /U. Thus, energy gain by 
•J'q for small n is ascribed only to the potential part even for small \U\/t. Anyway, 
energy is reduced from that of GWF for arbitrary values of n and U/t (U < 0). 

Shown in the insets of Figs. 15(a) and (b) are the optimized values of Et/t and 
Ejj/U respectively, the behavior of which is basically the same with that at half 
filling. We have estimated Uq for all the available densities of L = 10 in a similar 
manner to half filling, and showed it in Fig. 19. The value of Uq thus determined is 
weakly dependent on n: \Ug\/t increases subtly with decreasing n from half filling, 
has a maximum about 9.2 around n = 0.5, then slowly decreases and seems to tend to 
\U\ = W = 8t in the low density limit. This behavior is qualitatively similar to that 
of GA discussed in Appendix A. Incidentally, it is known that in the low density limit 
an s-wave onsite pair is formed for any negative value of U for the square lattice, 2 ) 
and is connected to s-wave superconducting states in finite densities. Nevertheless, 
accurate properties in this limit within the normal phase seem still unsettled. 

Fig. 19. Critical values of metal-insulator transition in 2D attractive Hubbard model for the up- and 
down-spin binding wave function and for Brinkman-Rice transition in Gutzwiller approximation. 
The magnitude of circle indicates possible errors. For the former system-size dependence is 
negligible in this scale; in fact the results for L = 12 (n = 0.5 and 1.0) are included among those 
for L — 10. The crosses indicate the critical value of metal-insulator transition by DMFT with 
a half-ellipse DOS. 8 ' The boundary U m of another instability arising in GA is also shown for 
2D by a dash-dotted line. In the low-density limit, Um/t ~ —4.9. Since chemical potential is 
constant with respect to n between the lines U m and £7br, the phase there is unstable. Since 
DOS of the square lattice is flat at the band edges, the region of this instability does not spread 
to the vicinity of half filling, in contrast to the case of a half-ellipse DOS. 37 ' The explanation 
as to Um and BRT will be given in Appendix A. 

Now, we go on to the correlation functions by &q. Shown in Fig. 20(a) is 
momentum distribution in both sides of Uq for some values of n. In every density 
definite discontinuities at kp can be seen in the weak correlation side (U/t = —8 
indicated by solid symbols), whereas in the strong coupling side (U/t = —12, open 
symbols) we can hardly find anomaly at k-p. Thus, the features near the Fermi 
surface are basically accords with those at half filling. In Fig. 20(b) N(q) and P(q) 
is depicted for some values of n. Here, although the data in the weak correlation 
side (|£7|/i = 8) is employed, ones for \U\ > \Uq\ make no difference qualitatively. 
As for N(q), as soon as electron density leaves half filling, N(G) greatly diminishes 
and the maximum position moves to an incommensurate wave number. Accordingly, 
dominant antiferro CDW correlation is restricted to half filling. In contrast, the peak 
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of P(q) at q = remains sharp even for n = 0.24; singlet-pairing correlation survives 
in the whole electron density. All these features are the same with those of GWF. 
Lastly concerning the gaps, the small- |g| behavior of N(q) and S(q) is essentially 
independent of n, thus the same with that at half filling. Accordingly, &q is always 
gapless for charge degree of freedom, whereas a finite spin gap opens in the insulating 
phase. 

Fig. 20. (a) Momentum distribution function by 'Z'q for four values of n and U/t — —8 and —12. 
(b) Structure factors of charge density and onsite singlet pairing by $q for four values of n and 
U/t = —8. Inset shows spin structure factor on the T-X line at half filling for the metallic (solid 
symbols) and insulating (open ones) regions, where the values of U/t are —9, —10, —11, —13, —16 
and —32 from above. In each panel 10 6 -2 x 10 6 samples are used for each parameter set. 

4.4. Discussion and comparison 

First of all, we compare the above results with those of DMFT. 7 )' 8 ) In DMFT 
there are two types of solutions — metallic (\U\ < |£7 C 2|) and insulating (\U\ > \U c i\) — 
in the normal phase, although details are different among different conditions or 
methods used to solve the Anderson impurity problem. In a intermediate regime 
l^cil < \U\ < \U C 2\, where two solutions coexist and a kind of phase separation may 
arise, a metal-insulator transition takes place. Indicated in Fig. 19 by crosses is 
that critical value, which seems to decreases monotonically with decreasing n. The 
fact that the values of VMC and DMFT broadly coincide corroborates the existence 
of a metal-insulator transition, which brings about a change in the mechanism of 
superconductivity, even in a realistic dimension. 

The above unstable phenomenon in the intermediate correlation strength had 
been reported for the half-filled RHM in a magnetic field by DMFT. 37 ) According 
to it, in the intermediate regime of U there appear a jump and a hysterises in the 
magnetization curve. This instability was considered as an analogue of the one 
arising in GA (equivalently GWF in this case), which we study in Appendix A. 
However, there is an evident difference in that the magnetization in DMFT does 
not saturate after this discontinuity. Here, we point out that this instability found 
in DMFT has nothing to do with the one in GWF. To begin with, recall that the 
range of instability in GWF necessarily includes n = 0, namely total energy is a 
linear function of n in a finite range near n = 0, as shown in Fig. 23. On the other 
hand, in Fig. 21 we show the same quantity by \Pq. In this case the total energy is 
a concave function of n in the low density regime for any coupling strength, namely 
1/Xc = d 2 E/dn 2 > 0. Thus, there is no instability for n ~ 0. It follows that the 
instability brought by GWF is spurious at least for small n. Incidentally, E/t for 
L = 10 shows a tendency to have a cusp-like winding at n = 0.6 (also slightly at 0.32) 
and be convex around n = 0.7 for L = 10, similarly for L = 8 a cusp at n = 0.5625 
and convexity at n ~ 0.7. This is probably caused by the finite system sizes, 38 ) and 
has nothing to do with the instability expected in the critical regime by DMFT. 39 ) 

Now, returning to Fig. 17, we compare Z of <Pq with that calculated from the 
self energy in DMFT. The behaviors of Z's roughly coincide at U — > and U ~ Uq, 
but does not for the intermediate values of U. It is probable that some intrinsic effect 
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Fig. 21. Energy of measured from the value at U = —oo (E x — Un/2) as a function of n for 
several values of U/t. Symbols denotes the VMC data, some of which are fit linearly in low 
density as shown by dashed lines. To check the system-size dependence in a high-density regime, 
the data for L — 8 are simultaneously plotted with those for L — 10. 

of electron correlation is not sufficiently introduced in for the metallic regime. 
Indeed, Z of &q is very close to that of GWF for small \U\. 

Then, we look at the improvement in E/t on GWF. Figure 22 shows the energy 
decrement by &q for some L's at half filling and simultaneously for some n's for 
L = 10. In every case, AE increases exponentially at first with increasing \U\/t, 
reaches a maximum around U/t = —12, then decreases as t/\U\. The improvement in 
the metallic regime is certainly small. The chief reason of the maximum at U/t ~ — 12 
is the abrupt change of E — Un/2 in GWF, as seen in Fig. 5. As for the transition of 
•Fq, the change of E — Un/2 around Uq is unexpectedly gentle and does not cause a 
conspicuous effect on AE. Moreover, the large system-size dependence is attributed 
mostly to GWF. Accordingly, one has to be deliberate in adopting E(&g) as a 
standard of the normal-state energy, for example, when one estimates a condensation 
energy. 

Fig. 22. Energy difference between GWF and \Pq as a function of \U\/t. At half filling data for four 
system sizes are plotted with solid symbols. The value of \E\/t itself decreases with increasing L 
for both $g and but the decrement of E(^>g) is thrice larger than that of E(^>q) in comparing 
L = 12 with L — 6. For U/t = —12, which gives the maximum, AE in the thermodynamic 
limit is extrapolated as 0.204 by second-order polynomial fit from the data up to L = 20 for 
GWF and up to L — 12 for <?q. Simultaneously, the data of three densities for n < 1 (L = 10) 
are shown with open symbols. The dashed [dash-dotted] line is an arbitrary one proportional 
to t/\U\ [exp(-t/U)]. 

§5. Summary 

In this paper we have studied normal-state properties of the attractive Hubbard 
model, especially in 2D, based on a series of variational Monte Carlo calculations. 
Main results are itemized as follows. 

(1) In the Gutzwiller wave function a precipitous quantitative change in various 
properties takes place around an interaction strength \U\ = \U co \ somewhat larger 
than the band width. In the region of\U\ > \U co \ the onsite pairing correlation grows, 
and the effective mass of electrons becomes large. In both sides of U co , however, a 
clear Fermi surface exists and there is no excitation gap. Hence, this change in 
GWF is a metal-to-metal crossover, and is distinguished from the Brinkman-Rice 
metal-insulator transition, although the value of U co is similar to Ubr- Furthermore, 
we have confirmed with better statistics that BRT is an artifact of the Gutzwiller 
approximation, and does not arise in finite lattice dimensions. On the other hand, 
we have mentioned GWF itself has unphysical features to be refined in correlation 
functions and the phase stability. 

(2) An improved wave function &q with a binding between j- and j-spins in 
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adjacent sites has been introduced by considering the perturbation expansion in 
the strong correlation limit. This wave function succeeds to the merit of GWF 
but improves its weak points; undergoes a substantial metal-insulator transition 
at \Uq\ roughly of the band width, which is favorably compared with the result 
of dynamical-mean- field approximations. We have constructed a phase diagram in 
the U/t-n plane; in the weak-coupling regime \U\ < \Uq\ Fermi surface is clearly 
recognized and there is no excitation gap, whereas for \U\ > \Uq\ the jump of n(k) 
at k = &f vanishes and there appears gap behavior in spin degree of freedom in 
accordance with QMC. 41 ) The state in the insulating regime is different from those 
of DMFT and of GA in that kinetic energy does not vanish. 

(3) We have considered a Jastrow-type long-range correlation factor as an al- 
ternative improvement. At half filling this wave function undergoes a first-order 
transition to the antiferro CDW at a finite value of U, whereas for incommensurate 
fillings, such a transition does not occur. In the metallic regime weak attractive in- 
tersite correlation as well as appropriate attractive onsite correlation slightly reduces 
the energy, but the state is essentially the same with GWF. Once a CDW state is 
realized, singlet pairing correlation is severely suppressed in this wave function. 

In the remainder, we make a couple of short discussions with future studies in 
mind. As for the metal-insulator transition the result of is consistent by and large 
with that of DMFT. Meanwhile, a recent study by a non-self-consistent T-matrix 
approximation 42 ^ concluded that in case of 2D a pseudogap in the single-particle 
spectral weight is induced even in a weak- interaction regime \U\ <C W due to the 
pair fluctuation, 43 -* if the temperature is sufficiently close to T c . In fact, a sign of 
pseudogap is also recognized for \U\ < W/2 in QMC data, 4 -*' 5 ) although this feature 
does not seem restricted to the 2D square lattice in seeing similar QMC results 
obtained for 2D triangular and 3D simple cubic lattices. 44 ) To study this issue &q 
has room for improvement, because pairing correlation as well as pairing fluctuation 
is not introduced deliberately for \U\ <C \Uq\ as mentioned in §4.3. 

Next, we refer to the difference between AHM and RHM. A common aspect 
among the three normal-state wave functions studied is that many properties in 
AHM only weakly depend on electron density except for the CDW correlation. This 
aspect is also common to a superconducting state in AHM, 46 ) but in sharp contrast 
with that of RHM, in which Mott transition is limited to the high-density regime 
especially at half filling. In addition, we have shown that !?j hardly improve for 
AHM unlike for RHM, and that for n < 1 the behavior of correlation functions by 
GWF is distinct between AHM and RHM. Moreover, in the insulating regime &q is 
gapless in charge degree of freedom, but is gapped in spin one for any electron density. 
This corresponds to the accepted conception that RHM at half filling has a Hubbard 
gap in charge part but is gapless in spin part. Nevertheless, if a similar insulating 
state is assumed for n < 1 in RHM, the gap behavior on spin and charge must be 
reversed from half filling. For this reason, it is not appropriate to discuss RHM in 
n < 1 on the analogy of AHM with h = 0; one should directly study RHM. Since <Pq 
is not properly extended to less-than-half filling for RHM, it is important to study 
a wave function in which intersite correlations are systematically introduced. 45 ) 

In this paper we have not touched on superconducting states, although in the 
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real ground state of AHM an s-wave superconducting order prevails all over the 
parameter plane, and an antiferro CDW order coexists particularly at half filling. 
Actually, trial wave functions with such orders further reduce the energy; we will 
treat this subject in a coming publication. 46 ) 
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Appendix A 

On Brinkman-Rice Transition 



In a previous paper 16 ) we discussed the absence of the Brinkman-Rice metal- 
insulator transition (BRT) for RHM with the VMC results. Afterward, the exact 
formula for GWF in D = 1 and D = oo (D: lattice dimension) were diagrammatically 
derived; 23 )' 47 ^ for D = 1 the absence of BRT is proven, whereas for D = oo the exact 
formula coincide with the GA formula, so that BRT exists. As for 2 < D < oo, there 
is no proof of the absence, but van Dongen et al. 47 ^ concluded the absence of BRT by 
making a couple of plausible assumptions about the scaling properties of spin-spin 
correlation function. In this Appendix we reconsider it with refined VMC data for 
U < 0, and refer to the dependence on lattice dimension. A part of the discussions 
on this transition was already mentioned in Ref. 48). 

First of all, let us look at the Gutzwiller approximation (GA) for U < 25 ) and 
H = briefly, which are obtained in a similar fashion to the repulsive case. 20 )' 21 ) 
According to GA, total energy is given by E = qe$ + Ud for n-t = n\ = n/2, where £q 
is the noninteracting energy per site for given n, and q is the renormalization factor 
in GA: 

d)(l-n + <£) + y/(%-d)d 



1 



. 2 



2 



(A-l) 



E has to be minimized with respect to d. With increasing \U\/t electrons more and 
more form onsite pairs, and singly occupied sites finally vanishes at a critical value 
U = Ubr, where d = n/2, q = and 7 = l/g = 0. Thus, electrons are completely 
localized and kinetic energy vanishes. In contrast to the repulsive case, BRT takes 
place at arbitrary electron density for U < 0. In Fig. 19 the critical value |J7br| 
for the square lattice is plotted together. At half filling the critical value has the 
same magnitude with the repulsive case: J7br,/£ = — 8|eo|/t (= — 12.969 •••); with 
decreasing n, | C/br | / 1 once increases slightly, has a maximum at n ~ 0.62, then 
decreases. The behavior for n — > is singular. 

Before examining BRT, we touch on another instability arising in GWF. In 
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connection with magnetic properties of 3 He and heavy-fermion systems, some papers 
treated RHM in a magnetic field H. Magnetization curve m(H) obtained by GA 21 ) 
has a discontinuous jump to the full moment at a value of H for U m < U < Ubr with 
U m = ££/br (£ = 0.38 for 2D square lattice). Namely, there exists an unstable region 
of m near m = 1 for intermediate values of U/t. This fact is interpreted for AHM 
through the canonical transformation eq. (2-2) as follows: For \U m \ < \U\ < \Uer\ 
there exists a finite region for small n where chemical potential £ is constant, namely 
particle number is variable. We show the boundary of this unstable region also in 
Fig. 19. Unlike BRT, this instability does not stem from GA but from GWF, because 
quantitatively similar jumps are observed in VMC results for H?q. 49 ) This instability 
is confirmed for AHM as the constant chemical potential [£ = (dE j 'dn)i] as shown in 
Fig. 23, where total energy of GWF by VMC calculations are plotted versus n. Total 
energy is a linear function of n for small n, and deviates roughly at GA's critical 
values indicated by arrows. As long as this instability is concerned, GA is broadly 
correct. Anyway, this instability near n = is an artifact of GWF and removed in 
<Pq as discussed in §4. 



Fig. 23. Energy of GWF measured from the value at U = —oo (Eoo = Un/2) as a function of n for 
several values of U/t. Symbols denote the VMC data, which are fit linearly in low density as 
shown by dashed lines. Arrows indicate GA's critical values of n (at Um) under which chemical 
potential is constant. 

Now, we return to BRT. Since this transition is caused by the behavior of kinetic 
and potential energies for small-7 region, we consider the expansion in 7, first. In 
GA the leading terms of both Et/t and Eu/U for 7^0 are linear in 7 for arbitrary 
n; at half filling, in particular, they are simply written as, 



t t(l + 7) 2 



= 4^7(1- 2 7 + ■■■). 



Eu 1 1 l n , , 
= = - 1 — 7 + • • •)• 

U 21 + 7 2 V ' ; 



The coefficient of the leading term of total energy E tot = E t + Ey changes its sign 
at I £7 1 = 8|eo|t. Thereby, it is determined whether the minimum is situated at 7 = 
(namely localized) or not. For the situation does not change. 

Now, let us consider GWF without approximations. For 7(= 1/g) — > 0, GWF 
is expanded as, 



*g = n [i - a - g)dj\ *f = ^ n b*i + d i\ ^ 

+ 7 ^e/2^iVe/2) 



g N 7 N-N c /2 



4 0) +7^+7 2 4 2) + 



(A-2) 
(A-3) 



with dj = rij^riji and dj 
sites* 

~(N-N c /2+i) (N e /2-e) 
3 k&j) 

where all the configurations with N c /2 — I d-sites are summed up. If we assume 
the expectation values of kinetic and potential energies by can be expanded with 



G 



1 — rij^riji. Here, $q is the state with 21 singly-occupied 

(A-4) 



E 

{d:N e /2-£} 
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respected to 7, they are formally written as 



(*gI*g) 



jKi + 7 2 K 2 + 0( 7 3 ) 



Po + 7 2 ^2 + 0(7 4 ) 



(A-5) 



(A-6) 



0Pg|#g) 

Here, ^ = Ti/JV , K 2 = T 2 /A r o, ^0 = ^0/^0 and P 2 = (C/ 2 - UqN^Nq)/^, with 



(vW\Ht\Vg>) + (Vg'lHtl^) /t, T 2 



(^ /2) \Hu\V^)/U, 



Henceforth, we concentrate on the half-filled band, because system-size dependence 
can be checked and the situation does not change for n < 1, in addition to another 
merit that the discussion here is directly connected to the conventional repulsive case 
with n = 1 and h = through eq. (3-4) 

For finite systems the situation is simple, because the expansion coefficients K^s 
and P/s are finite, namely the expansion by 7 eqs. (A-5) and (A-6) are pertinent. 
This is seen in Figs. 24(a) and (b), where E t /t and Ejj/U are plotted versus 7 (~ 0) 
for some values of L. Since the leading terms of E t /t and Ejj /U for 7 — ► are linear 
and quadratic, respectively, E tot has a minimum at 7 = — tifi/(C/P 2 ) for \U\ 3> t. 
Thus, 7 does not vanish as long as \U\/t is finite: BRT is absent. 



Fig. 24. (a) Kinetic and (b) potential energy of GWF as a function of 7 near 7 = 0. System-size 
dependence is compared. The results of GA are also shown by a dash-dotted line. In (a) the 
leading linear term for L = 00 is shown by a shadowed line. Its slope K\ (See Table I) is 
probably a half of the value by GA: 4eo/t. The dashed line smoothly connects the data for 
L — 20. Error bars are attached only for L — 20 and 24 for clarity. In (b) the data for L = 6-20 
are fitted by quadratic curves (solid line). As L increases, the range of 7 where Eu /U obeys 
n/2 — const, x j 2 shrinks; dashed lines connecting data points for L = 24 are guide for eyes. 
5 x 10 4 -10 7 samples are used. 

Nevertheless, the situation is not so simple for the thermodynamic limit. Note 
that in Fig. 24 the behavior of Et/t for 7 ~ seems almost independent of L, whereas 
that of Ejj/U is highly dependent. Therefore, let us consider the coefficients K±, K2 
and P 2 , which are also obtainable by VMC calculations. In Fig. 25(a) the coefficients 
thus obtained, which are naturally independent of 7, are shown for L = 10 and for 
some values of 7. The values thereby and similarly for different L's are plotted as a 
function of 1/L in Fig. 25(b). K\ only weakly depends on L and is well fitted linearly, 
whereas K2 and P 2 largely depend on L and diverge as —a — b/{l/L) c for 7^0. 
It follows the expansion form eqs. (A-5) and (A-6) is inadequate. This situation by 
VMC analysis is quite common to the ID case. 4S ) On the other hand, according to 
the exact analytic results for the ID repulsive model, 23 ^ the leading term for Eu/U 
has logarithmic correction as g 2 ln(l/g), while Et/t = 2geo/t for n = 1 and g — > 0. 
Thus, it is probable that such a nonanalytic factor appears also in higher dimensions. 
Of course, logarithmic correction does not affect the absence of BRT, though. 
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Fig. 25. (a) Low-order expansion coefficients of 2D GWF for 7 — > and n = 1.0 calculated by 
a VMC method. Incidentally, P = n/2. The least square fit yields = -3.24 ± 0.01, 
K2 = —30.90 ± 0.09, and P2 = —5.83 ± 0.02. Accuracy of calculations becomes high around 
the maximum of l^o'l 2 , e.g. 7 ~ 0.045 for L — 10 as seen in the inset of (b). (b) System-size 
dependence of low-order expansion coefficients. K\ for L = 00 is estimated as —3.235 by the 
second-order polynomial fit. The other curves are well fitted by the formulae —a — b/(l/L) c , 
where a = 3.79, b = 0.932, c = 1.46, for K 2 and a = 0.86, b = 0.173, c = 1.46 for P 2 . The power 
c seems common to K2 and P2 irrespective of lattice dimension (see Table I). Inset shows the 
weight of <p£\ that is 7 2 (*-^)(#£Vg V^gI^g), as a function of 7. 

Instead of the expansion, let us consider Eu/U directly, which is not a divergent 
quantity. Assume the leading power of Eu/U is continuous: Eu/U=E oc /U— const, x 
7 a . Then, the optimized value of 7 is nonzero for finite U/tiia > 1, since the leading 
term of Et/t is linear in 7 due to the above discussion. In Figs. 26(a) and (b), Et/t 
and Eu/U are depicted for several small values of 7 as a function of 1/L 2 . Since the 
values changes rapidly with increasing L and we do not know a priori the system-size 
dependence, reliable extrapolation is not anticipated from the raw data. 50 ) Therefore, 
we consider the ratio between the values of different 7; here we use 7 = 1/201 as a 
reference. In Figs. 27(a) and (b) such ratios are plotted for E t /t and (E^ — Eu)/U, 
respectively. In this case Eqo/U = 0.5. In this plot the system-size dependence is 
smaller, so that the behavior in the thermodynamic limit can be estimated to some 
extent. The crosses on the vertical axis indicate the values if Et/t or (-E00 — Eu)/U 
is precisely proportional to 7, namely 1/(2017). Similarly, the dots indicate the case 
of 7 2 : l/(20l7) 2 . In Fig. 27(a) all the curves substantially point to the crosses; E t /t 
is proportional to 7, as expected. On the other hand, in Fig. 27(b) the extrapolated 
values seems by far nearer to the dots than to the crosses especially for small values 
of 7; the effective power a ~ 2 (at least > 1). Some deviation from purely linear 
or quadratic behavior for large L may be attributed to possible nonanalytic factors, 
besides using small but finite magnitude of 7. Thus, we can confidently conclude 
BRT is absent not only in ID but in 2D. 

Fig. 26. System-size dependence of (a) E t /t and (b) Eu/U for several small values of 7 on the 
square lattice. Symbols are common in (a) and (b). 5 x 10 5 (L = 24)-10 7 (L = 6) samples are 
averaged. 



Fig. 27. System-size dependence of (a) E t /t and (b) [-El/ (7 = 0) — Eu]/U for several small values of 
7 normalized by the value at 7 = 1/201. Crosses (dots) on the vertical axis indicate the values 
assuming E t (-y)/t or [E^ — Eu(j)]/U is precisely proportional to 7 (j 2 ). 

Similar analyses have been carried out for the simple cubic lattice. Compared 
with 2D in Fig. 26, the behavior of E t /t and Eu/U in 3D (not shown) is qualitatively 
similar, but the rapid change for 7^0 occurs at a smaller value of L (or a large 
value of 7) in both energy components. In Fig. 28 the ratios corresponding to Fig. 27 
are shown for 3D. Although the feature is substantially the same with that in 2D, 
smaller 7 is needed to single out the behavior of leading terms. To confirm this 
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aspect quantitatively, we compare the expansion coefficients among ID, 2D and 3D 
lattices in Table I. As for the kinetic part, K\ scarcely changes with the system 
size in each dimension, and the relation K\ = 2eo/t holds not only for ID, but 
probably for other finite dimensions (so/t = — 4/7T, — 16/-7T 2 , —2.0048 • • • for D = 1, 2 
and 3, respectively). This point is in sharp contrast to the 1/D expansion discussed 
below. On the other hand, with increasing lattice dimension not only the absolute 
values of P2 and K2 increase, but their divergent behavior for L — > 00 becomes more 
remarkable. For instance, compare the power c in Table I. Consequently, each energy 
component by GWF quantitatively approaches the value of GA with increasing D, 
but the qualitative features for 7 (g) — > remain distinct. 

Fig. 28. System-size dependence of (a) E t /t and (b) [Eu('"f = 0) — Eu]/U normalized by each value 
at 7 = 1/201 for several small values of 7. Crosses and dots on the vertical axis indicate the 
same with Fig. 27. Samples as many as 2.5 x 10 5 (L = 10)-5 x 10 6 (L = 4) are averaged. 

Table I. Expansion coefficients of Et/t and Eu/U for ID, 2D and 3D lattices. In the last line but 
one the values of 2eo/t are entered for comparison with K\. Entered in the last line are the 
common powers c for P2 and K2 in the fitting form of —a — bL c . For ID, the data in ref. 48 ' for 
RHM are adopted; the values are adjusted to the attractive case through eq. (3.5): Ki = Kg and 
P 2 = — P2 from the original data (Ke, P 2 ). They are estimated from fewer samples (5 x 10 4 -10 5 ) 
with periodic boundary condition. 



ID chain 


2D square 




3D simple 


cubic 
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0.12 
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1.54 
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2.76 



This point has been already argued by the 1/D expansion analysis. 47 ^ In this 
approach the starting point is D = 00, where GA becomes the exact estimation of 
GWF, and properties for finite dimensions are improved by 1/D (and 1/D 2 etc.) 
corrections on GA. Owing to this correction, most properties are quantitatively im- 
proved even in ID. However, the exception is the behavior for 7^0. In 1/D 
expansion, the leading terms of the energy components for g — > are calculated as 
Et/t = c\gso/t and Ejj/U = C2g, where 




This indicates that 1/D corrections only reduce the coefficient and Ejj/U remains 
linear in g for small g. Therefore, this approach does not get rid of BRT, but only 
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increases the value of |E/br|- In this connection it was pointed out that D = oo 
is singular as to the properties for g — > 0, so that 1/D expansion is inappropriate 
as long as near BRT. 47 ) Similarly, the coefficient c\ is reduced by 1/D corrections, 
which means different values of c\ for different -D's. However, the value of C\ by 
VMC is constant, that is 2, for the three realistic dimensions. Thus, c\ also seems 
singular at D = oo. 

Summarizing, BRT is restricted to the infinite dimension, and never takes place 
in finite dimensions. D = oo is singular for the properties in 7 (g) — > 0; the linear 
coefficients in the 7-expansion of both E t /t and Ejj/U have discontinuity at D = 00, 
specifically as 2eo/t — ► 4eo/t and — > 1/2, respectively. It is a remaining issue how 
logarithmic corrections, which are not relevant in ID, affect the power counting, if 
any for 2 < D < 00. 

Appendix B 

Possibility of long-range orders in GWF 

In this Appendix we study the singular behavior of correlation functions calcu- 
lated with GWF as a supplement to §3.3. To check the system-size dependence we 
have plotted the inverse of N(G), P(0) and S(G) versus 1/L 2 (L = 6-20) for sev- 
eral values of 7 (not shown) , and found that they are smoothly extrapolated toward 
L = 00 by polynomials. Figure 29 shows the values thereof as well as the original 
finite-size data versus 7. From this figure we can read a couple of facts. First, both 
N(G) and S(G) have very little system-size dependence for large values of 7, while 
they come to have appreciable size dependence with decreasing 7. The value of 7 
at which system-size dependence becomes remarkable approximately coincides with 
U = U co (7 ~ 1.2); simultaneously the system-size dependence changes its tendency 
from convergent to divergent (not shown). Thus, the crossover of the states also 
affects the correlation functions. Second, although N(G) or equivalently P(0) does 
not diverge for finite coupling strength, in the 7^0 limit the extrapolated inverse 
quantity of L — > 00 is likely to vanish. Namely, N(G) and P(0) diverge in the limit 
of 7 = or \U\/t = 00. On the contrary, the order parameter Ocdw defined in 
§2.1 measured in the same VMC sweep is always zero within statistical errors for 
any values of 7 and L. From these results we reach a conclusion that GWF does 
not possess CDW and singlet superconducting long-range order even in this limit. 
This point resembles the feature of ID Heisenberg model. Finally, concerning S(q), 
it is natural that the structure in S(q) vanishes and its amplitude fades out with 
increasing \U\/t, according to the diminution of isolated spins. 

Fig. 29. Dependence of (a) 1/N(G) [or equivalently 1/2P(0)] and (b) S(G) on a variational pa- 
rameter 7 at half filling for various system sizes and extrapolated values for L = 00, which are 
obtained by polynomial fit. Insets are magnification for small-7 region. Solid lines are fitted 
curves by proper polynomials. Dashed lines for L — 00 are guide for eyes. 2.5 x 10° (L = 20)- 
I0 6 (L = 6) samples are used for each parameter set. 

To see the system-size dependence in less-than-half filling, we take up a quarter 
filling. As seen in Fig. 10(a), N(q) is enhanced for any q by attractive 7 also in this 
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density. On the path r — > X — > M — > T there are three singular g-points (A, B and 
C), which correspond to 2k? or 2{G — 2&p) in the noninteracting system indicated 
by arrows. In contrast to half filling, however, the increment at these singular points 
is fairly small. 51 ) We plot 1/N(q) at these points versus 7 for three system sizes 
in Fig. 30. Although system-size dependence is not simple, it is certain that N(q) 
converges even at 7 = for each of the three singular wave numbers. Thus, the CDW 
correlation in GWF is not intrinsic even at a quarter filling, where the band filling is 
quasi commensurate. On the other hand, P(q) in Fig. 10(b) is still greatly enhanced 
at the r point by attractive correlation. As shown in Fig. 30, the behavior of P(0) for 
less-than-half filling is essentially the same with that for n = 1 [Fig. 29(a)]; although 
P(0) diverges in the limit of 7 — > and L — > 00, a long-range order of singlet pairing 
is not realized. The value of U at which appreciable system-size dependence appears 
is approximately U co (7 ~ 0.1). Consequently, the behavior of N(q) and P(q) of 
GWF basically matches up to the accepted phase diagram of AHM. 

Fig. 30. Inverse of N(q) at three singular wave numbers q = Qa, qB and qc and inverse of P(q) at 
q = (0,0) as a function of 7 for three system sizes at a quarter filling. Concerning P(q), the 
values at half filling are also shown for comparison. The singular wave numbers of N(q) are 
incommensurate and subtly different size by size; thereby the system-size dependence of N(q) 
is more complicated. 



Appendix C 

Jastrow-type long-range correlation factor 

In this Appendix as an alternative prescription to improve the normal state we 
introduce a distance-dependent correlation factor. This orthodox Jastrow-type wave 
function was useful for RHM in n < l. 27 ) We describe the wave function in §C.l, 
and discuss the VMC results in §C.2. 

C.l. Jastrow-type correlation factor 

It is natural that even in the model with only onsite interaction the influence 
of interaction is not confined to the same site. In RHM the importance of repulsive 
intersite correlation factors were obvious, because configurations in which electrons 
are mutually apart are favorable for both kinetic and potential terms. Nevertheless, 
in AHM the situation is not so simple, because electrons in proximity seems, if 
anything, advantageous to the potential term but at a glance disadvantageous to 
the kinetic term. The aim of this section is to study how a simple long-range factor 
works in such a system. 

In this work we adopt a two-body Jastrow-type wave function with a spin- 
independent correlation factor: 

Vj = Vj<Pf= II {l-[l-r)(rij)}nian jT }$ F , (C-1) 

i>j,UT 
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where iy = rj — rj. For rj(r) we employ one of the simplest forms: 

( g (r = 0) 

" (r) = \ [v^K^^ (r^O) ' (C ' 2) 

with r = (x,y). Such a form was originally introduced to describe a Tomonaga- 
Luttinger liquid for RHM and the t-J model in ID, 52 ) and later extended to 2D. 
53 ) Nevertheless, we overlook the subtle difference in long-range behavior between 
Tomonaga-Luttinger and Fermi liquids, 54 ) because our concern here is how a distance- 
dependent correlation factor works. The wave function $j with eq. (C-2) has two 
variational parameters g and v. The former one g or 7 (= 1/g) controls double occu- 
pancy as that in Vq, and encourages (discourages) the increase of d for 7 < 1 (7 > 1). 
On the other hand, r/(r) substantially works as r v in a short distance (r <^ L), so 
that i](r) plays a repulsive or an attractive role for the intersite part according as 
v > or v < 0. When v is zero, the intersite part in eq. (C-2) becomes constant and 
$j is reduced to GWF. Incidentally, $j for U < is not simply connected to the 
counterpart for U > through eq.(2-2) even at half filling 55 ) unlike GWF and )Pq. 

An additional interest in <Pj is whether a phase separation or a CDW arises 
in AHM. For, when this type of wave functions is applied to the t-J model, 52 ) " 54 ) 
which has an attractive parameter J, phase separation is brought about as a first- 
order transition at a finite value of attractive interaction, which is quantitatively 
accurate as compared to other reliable results. 

C.2. Results of VMC calculations 

Now we start with the half-filled band. In Fig. 31 expectation values of the energy 
components are shown for a variety of parameter sets. At first let us look at the case 
of attractive onsite correlation (7 < 1). As for the kinetic part, \E t \ is suppressed for 
strong attractive values of v ( ^ — 0.3) as anticipated and rather unexpectedly for 
v > on account of the tendency toward a fixed configuration with the maximum 
interparticle distance. Nevertheless, a weak attractive intersite factor (—0.25 ^ v < 
0) is somewhat advantageous to the kinetic term. This trend stems from the hopping 
between double occupied sites and neighboring empty ones induced by the introduced 
attractive intersite factor. This is supported by the behavior of d (= Ejj/U), which 
is slightly reduced for the corresponding small negative values of v. On the other 
hand, for large negative values of v, d is enhanced as expected. Furthermore, note 
that repulsive intersite correlation also promotes double occupation. 



Fig. 31. Expectation values of (a) kinetic and (b) potential energies of tf'j at half filling as a function 
of onsite correlation parameter 7 for some values of intersite correlation parameter v. In (b) 
data for v — —0.1 and v = —0.25 almost overlap. To survey a wide range of parameters, data 
of a small system L = 6 are shown; there is no qualitative difference from larger systems. 10 6 - 
2 x 10 6 samples are averaged for each parameter set. Statistical errors are much smaller than 
the symbol size. 



In Fig. 32 total energy E/t is plotted for three values of U/t. Figure 32(a) 
shows a typical case of small \U\/t, where the energy minimum is given by (7, v) = 
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(0.524,-0.05). This optimal state is realized by small attractive intersite corre- 
lation with a compensatory enhancement in the attractive onsite part. Since the 
improvement in energy on GWF (0.548,0) is very small as seen in the inset, this 
state is no more than a weakly perturbed GWF. In this case E/t monotonically 
increases with z/s going away from the optimal value —0.05. Figure 32(b) shows the 
case for U/t = —6.8, in which the minimum is given by a similarly small negative 
value of v (0.314, —0.05) indicated by 'Metal', therefore the realized state is still 
a weakly perturbed GWF. Nevertheless, a distinct feature from (a) shows up; E/t 
becomes non-monotonic with respect to v and there appears another local minimum 
at (0.851,0.75) indicated by 'CDW. These parameter values are quite different from 
those of 'Metal', namely weak onsite attraction and fairly strong intersite repul- 
sion. Depicted in 32(c) is the total energy for a slightly stronger value U/t = —7, 
where the optimal parameter set is switched to 'CDW (0.952,0.83) from 'Metal' at 
(0.299,-0.05). 

Fig. 32. Expectation values of total energy at half filling for three values of U /t, (a) —4, (b) —6.8 
and (c) — 7. Inset in (a) is the magnification of the minimum part. Arrows indicate minima, 
particularly in (b) and (c) large (small) arrows means the global (local but not global) minima. 
'Metal' and 'CDW indicate the attributes of the minima. Symbols are common to all the panels. 

In Fig. 33 we show the optimized parameters for > U/t > —32 and for three 
values of L. For \U\/t ^ 6, only onsite attractive correlation grows with v keeping 
almost a constant value —0.05. At U/t = Ucuw/t ~ —6, however, it switches to a 
distinct curve, along which both correlation factors become more repulsive with in- 
creasing \U\/t, and sizable system-size dependence appears. Thus, some first-order 
transition is anticipated at this critical value. Recall that the intersite repulsive 
factor enhances double occupation [Fig. 31(b)], but obviously hates a gathering con- 
figuration like phase separation. In fact, a transition arising here is from metallic to 
("7r, 7r)-CDW, as will be identified below. It is interesting that in the model with only 
attractive interaction the state can be stabilized by the competition among repulsive 
correlations. 

Fig. 33. Optimal parameter values of tf'j at half filling for three system sizes. Due to discrete 
parameter values used, the data sometimes show zigzags. Some of the corresponding values of 
U/t is shown. The value of \Ucr>w\/t decreases a little with increasing L. Inset shows energy 
improvement by ^Pj on GWF as a function of U jt at half filling for three system sizes. In the 
thermodynamic limit the maximum value of AE/t at half filling (at U/t ~ —12) is roughly 
extrapolated as 0.26. For \U\ > \U co \, AE/t is approximately proportional to t/\U\. 

Before discussing the correlation function, we add a couple of notes on energy. 
First, by studying the optimized energy components (not shown), one finds that 
both Et/t and Ey jU have such jumps at U = Uqdw as \E t /t\ decreases but \Ejj/U\ 
increases, in accordance with ordinary phase transitions to some ordered states. 
Second, as depicted in the inset of Fig. 33 the energy gain by $j remains very little 
for \U\ < \Uqdw\ [Fig. 32(a)], whereas AE abruptly increases for \U\ > |C/cdw|- 
This contrasts with a mean-field-type CDW wave function, 46 ) in which a transition 
takes place at U = 0, and thence AE is smooth. Third, the maximum of AE is 
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situated at U ~ U co of GWF, and the system-size dependence of AE is conspicuous 
for | U | > |C/cdw|- They are mostly attributed to GWF, as we have experienced for 

Now, we proceed to the correlation functions. As shown in Fig. 34(a), momentum 
distribution n(k) [eq. (3-6)] for \U\ < I^cdwI represents a typical metallic state with 
discontinuity at kp. In this wave function an unphysical feature of GWF described 
in Appendix B is remedied by a small but effective intersite correlation factor jjl = 
—0.05. On the other hand, since the discontinuity at kp almost disappears for 
\U\ > \Ucdw\, the state becomes insulating. Next, from S(q) in Fig. 34(b) it is 
found that spin correlation is suppressed especially for \U\ > \Uqdw\, which reflects 
the sudden increase of onsite spin pairs. In the limit of \q\ — > 0, S(q) seems a linear 
function of \q\ even for \U\ > \Ucuw\- Therefore, $j does not have a gap in the spin 
degree of freedom. Concerning onsite singlet-pair correlation [Fig. 34(d)], the sharp 
peak at q = 0, which is enhanced with increasing \U\ for \U\ < \Uqdw\, is suddenly 
suppressed at U = UqdWi an d P(q) approaches a constant 0.5 for \U\ > \Uqdw\- 
This indicates that for \U\ > \Uqbw\ pairing correlation is restricted to a very short 
range, especially onsite. 

N(q) [Fig. 34(c)] increases as a whole with increasing \ U\ in the weak correlation 
phase (\U\ < |C/cdw|)) whereas it decreases for \U\ > |J7cdw| except for q = G; N(G) 
becomes overwhelmingly large in the strong-coupling phase, suggesting an antiferro 
CDW order in this phase. To confirm the formation of a long-range order, we 
plot N(G) versus \U\/t for a couple of system sizes in the inset of Fig. 34(c). For 
\U\ > |J7cdw|) N(g) is proportional to the system size, and approaches the value 
of completely polarized CDW state, that is N, as \U\ increases. In fact, antiferro 
CDW patterns are evidently seen in all the snapshots taken during Monte Carlo 
sweeps in the strong correlation phase. Thus, we have positive proof of a CDW 
order. Incidentally, one observes in Fig. 34(c) that N(q) in the CDW phase is not 
a linear function of q around the F point but qP((3 > 1). Thus, a gap opens in the 
charge sector. 

Fig. 34. (a) Momentum distribution and correlation functions of (b) spin, (c) charge density and 
(d) onsite singlet pairing, calculated with at half filling. Path in the momentum space is 
basically the same with Figs. 7(a) and 8. For this system (L = 10) C/cdw ~ —6.3. The data for 
\U\ < |(7cdw| are shown by solid symbols, otherwise by open symbols, common in all the panels. 
For each value of U/t, 5 x 10 5 samples are used. Inset of (c) shows N(q) at q = G = (tv,tt) as 
a function of coupling strength for three system sizes. Crosses with arrows on the right vertical 
axis indicate the values for the completely polarized CDW: N. In addition, the order parameter 
Ocdw for L = 10 is drawn by a dashed line. 56 - ) 

Finally, we look at less-than-half filling, where antiferro CDW correlation loses 
its advantage, because the nesting condition is no longer satisfied. As for E t and 
Eu, the behavior is analogous to that at half filling (Fig.31); weak attractive intersite 
correlation and proper attractive onsite correlation slightly stabilize the energy. A 
clear difference from half filling is that such a Fermi-liquid state is stable even for 
very large values of \U\/t at least up to 32 and a CDW phase does not appear. 
In Figs. 35(a) and (b) energy expectation values for U/t = —16 are depicted for 
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n = 0.72 and 0.24, respectively. In contrast to Fig. 32(c) the minimum of 'CDW 
never shows up in both densities. 57 ^ Thus, we conclude that \Pj remains a metallic 
state for all the values of U/t, but improves GWF only a little. 

Fig. 35. Expectation values of total energy in partially filling (a) n = 0.72 and (b) n = 0.24 for a 
relatively large value of \U\/t. Arrows indicate energy minima. 5 x 10 5 -10 6 samples are used for 
each parameter set. 

As a metallic state the improvement by an orthodox Jastrow-type wave function 
is considerably small for AHM regardless of electron density. This results sharply 
contrasts with those for RHM, for which extensive improvement is achieved for n < 1. 
This is mainly because the distance dependence of electron correlation in AHM is 
not monotonic like the correlation factor in <Pj. An extension to a many-parameter 
algorithm 58 ^ is desired. 

Appendix D 

Expansion of &Q near Uq 

For a further discussion on the metal- insulator transition in <Pq, it is useful to 
consider the expansion in 1 — fi again [see Appendix C of ref. 27)]. In the vicinity of 
[X = 1, one can expand <Pq with respect to 1 — fi(= m) as, 

* Q = «4 0) + mV^ + m 2 ^ + -.., (D-l) 

*Q = [IK 1 " Q$ *G, *Q = {EQil IK 1 " Qi)\} *G, 

* 3 i&j) 

*q={ E QiQki II (i-O0]}afe. 

Thereby, the energy expectation value is expanded, by abbreviating \^ < q) to as 

e]® = (o\h\o)/n , e q 1] = [(o\n t \i) + (i|Wt|o)] /n , 

E Q 2) = [(2|Wt|0> + {l\H t \l) + (0|Wt|2) + (l\Hu\l) + M(0|W|0>] /JV , (D-3) 
with = (£\£). According to eq. (D-2), if Eq has a minimum at fi < 1, namely if $q 
is metallic, £7q ^ must be negative and finite. Actually, we numerically confirmed that 

Eq ^ is negative and its absolute value is increasing with L for ID systems; 27 ) this 
is probably true in the present 2D case. Thus, there is no metal-insulator transition 
within \Pq at finite U in the exact sense. 

Leaving aside this fact, we discuss the relevance of the leading order \Pq^ to the 
transition, consists only of such configurations as ^4 = (• • • Q T I t OIO'")i 
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in which at least one <r-spin site is adjacent to more than one — cr-spin sites (O 
indicates either d-site or e-site). This configuration is included in a higher-order 
contribution in 7 expansion than ones belonging to ^ like B = (• • • O O O T O I 
O " " ■)> so tnat wrtn decreasing 7 the weight of <Pq^ becomes smaller by order of 7 
than that of ■ Thus, the absolute value of E^q may become negligibly small as 
compared with Eq\ to which lower-order terms in 7 contribute [the last term in 

eq.(D-3)]. Actually, a conspicuous reduction of \Eq^\ with decreasing 7 is seen for a 
ID system in Fig. 20 of ref. 27). Thus, the leading term in Et/t becomes practically 
quadratic in the fx — ► 1 limit for relatively small values of 7 [e.g. 7 = 0.25 in the 
inset of Fig. 12(a)], in contrast to the linear behavior for large 7. Recall that the 

transition in <Pq takes place at 7 ~ 0.2 (Fig. 14). Consequently, the problem of 

(2) 

our concern virtually rests on Eq . Actually, as shown in Fig. 13, total energy E/t 
for large \U\/t does not look linear in 1 — /j, near /j, = 1 but becomes very flat as a 
function of /x, and the optimized energy given by \i < 1 is almost the same with the 
insulating case \i = 1. 

In conclusion, since the contribution of E^ is irrelevant at least for small values 
of 7, a substantial metal-insulator transition can take place even in <Pq. Nevertheless, 
a metal-insulator transition will be treated more clearly, if we further improve the 

wave function by allowing for the following point. Namely, in the standpoint of t/U 

(£) 

expansion all the configurations in with odd t should be classified in the higher- 
order terms than £-th. For example, the above configuration A (I = 1) should be 
assigned to the same order with the one B (1 = 2). 
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